R for Chemistry Data Analysis and Chemometrics
Visual Analytics in R
Jordi Cuadros, Vanessa Serrano
January 2022

Package Management

Installation and Loading of R Packages

R has many contributed to solve specific needs and problems in several disciplines. These are usually hosted in repositories, such as CRAN (https://cran.r-project.org/), Bioconductor (https://www.bioconductor.org/), R-Forge (https://r-forge.r-project.org/), ROpenSci (https://ropensci.org/), Stan packages (http://mc-stan.org/r-packages/)… Yet other packages are hosted in version control websites as GitHub or GitLab.

In all cases, to make use of a package, it has to be installed and loaded in memory.

installed.packages()[,1] # Lists the installed packages
(.packages())            # Lists the packages currently loaded

Some relevant resources to find useful packages are https://www.rdocumentation.org/ and https://cran.rstudio.com/web/views/.

To install a package, we use

install.packages("dplyr")

To load a package in memory, we use

library("dplyr")

Alternatively, RStudio offers a menu-based package management feature.

In a script, to avoid repeated installations of a package and to ensure its availability, we can use

if(!require("dplyr")) {
  install.packages("dplyr",
                   repos="https://cloud.r-project.org/")
  library("dplyr")
}

If the package is installed but not loaded, we can still access its functions by using the ::operator.

dplyr::band_instruments

To browse the help files for a packages, run

help(package="dplyr")

Additionnaly some packages have a vignette and some demo code that can be inspected with thevignette() and demo() functions.

Last, to learn more about packages and package development, you may want to read https://r-pkgs.org/index.html.

Data Import and Export

Data Files

It is often useful to save and load data for its transfer among different application or just to keep a copy of it for back-up or later use.

Although can read and write many different types of files, we will focus on the most relevant format for statistics, chemistry and chemometrics. We will discuss here how to handle

  • text files,
  • structured data files, such as XML and JSON,
  • files from other computational programs (i.e. Excel and Matlab),
  • RData files, the default format for R, and
  • some data files which are specific for chemistry and chemical analysis.

Text Files

Let’s start from the wine data set in the FactoMineR package.

if(!require("FactoMineR")) {
  install.packages("FactoMineR",
                   repos="https://cloud.r-project.org/")
  library("FactoMineR")
}
data("wine")
help("wine")
str(wine)
## 'data.frame':    21 obs. of  31 variables:
##  $ Label                        : Factor w/ 3 levels "Saumur","Bourgueuil",..: 1 1 2 3 1 2 2 1 3 1 ...
##  $ Soil                         : Factor w/ 4 levels "Reference","Env1",..: 2 2 2 3 1 1 1 2 2 3 ...
##  $ Odor.Intensity.before.shaking: num  3.07 2.96 2.86 2.81 3.61 ...
##  $ Aroma.quality.before.shaking : num  3 2.82 2.93 2.59 3.43 ...
##  $ Fruity.before.shaking        : num  2.71 2.38 2.56 2.42 3.15 ...
##  $ Flower.before.shaking        : num  2.28 2.28 1.96 1.91 2.15 ...
##  $ Spice.before.shaking         : num  1.96 1.68 2.08 2.16 2.04 ...
##  $ Visual.intensity             : num  4.32 3.22 3.54 2.89 4.39 ...
##  $ Nuance                       : num  4 3 3.39 2.79 4.04 ...
##  $ Surface.feeling              : num  3.27 2.81 3 2.54 3.38 ...
##  $ Odor.Intensity               : num  3.41 3.37 3.25 3.16 3.54 ...
##  $ Quality.of.odour             : num  3.31 3 2.93 2.88 3.36 ...
##  $ Fruity                       : num  2.88 2.56 2.77 2.39 3.16 ...
##  $ Flower                       : num  2.32 2.44 2.19 2.08 2.23 ...
##  $ Spice                        : num  1.84 1.74 2.25 2.17 2.15 ...
##  $ Plante                       : num  2 2 1.75 2.3 1.76 ...
##  $ Phenolic                     : num  1.65 1.38 1.25 1.48 1.6 ...
##  $ Aroma.intensity              : num  3.26 2.96 3.08 2.54 3.62 ...
##  $ Aroma.persistency            : num  2.96 2.81 2.8 2.58 3.3 ...
##  $ Aroma.quality                : num  3.2 2.93 3.08 2.48 3.46 ...
##  $ Attack.intensity             : num  2.96 3.04 3.22 2.7 3.46 ...
##  $ Acidity                      : num  2.11 2.11 2.18 3.18 2.57 ...
##  $ Astringency                  : num  2.43 2.18 2.25 2.18 2.54 ...
##  $ Alcohol                      : num  2.5 2.65 2.64 2.5 2.79 ...
##  $ Balance                      : num  3.25 2.93 3.32 2.33 3.46 ...
##  $ Smooth                       : num  2.73 2.5 2.68 1.68 3.04 ...
##  $ Bitterness                   : num  1.93 1.93 2 1.96 2.07 ...
##  $ Intensity                    : num  2.86 2.89 3.07 2.46 3.64 ...
##  $ Harmony                      : num  3.14 2.96 3.14 2.04 3.64 ...
##  $ Overall.quality              : num  3.39 3.21 3.54 2.46 3.74 ...
##  $ Typical                      : num  3.25 3.04 3.18 2.25 3.44 ...

To write a delimited text file from a data frame, we use the write.table function or any of its derivative functions. To read, read.table is to be used.

write.table(wine, file="wine.csv", sep=",", dec=".",
            quote=TRUE, fileEncoding="UTF-8", row.names=FALSE)
wine2 <- read.table("wine.csv", sep=",", dec=".",
           quote="\"", fileEncoding="UTF-8", header=TRUE)

In case we want a tab-delimited text file, we will use \t for the argument sep.

Other options exist for managing text-files.

In case the format is not well-defined, we may want to use readLines and writeLines to read or write the file line by line.

If the file is large, read_table from the readr package or vroom from the vroom package are usually faster.

Structured Text Files – XML

An example of an XML file containing the structure of caffeine can be found at https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/2519/record/XML. It can be read in R with read_xml from the xml2 package.

if(!require("xml2")) {
  install.packages("xml2",
                   repos="https://cloud.r-project.org/")
  library("xml2")
}

Information can be extracted using XPath selectors, https://www.w3schools.com/xml/xpath_intro.asp.

xmlObj <- read_xml(
  "https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/2519/record/XML",
  options = "RECOVER")
xml_ns_strip(xmlObj)  # Important!
xmlObj_Name <- xml_find_all(xmlObj,xpath=
  ".//PC-InfoData//*[text()='Traditional']/../../..//PC-InfoData_value_sval/text()")
as.character(xmlObj_Name)
## [1] "caffeine"

NMR spectra can be obtained in CML (Chemistry Markup Language, http://www.xml-cml.org/spec/) format from https://www.nmrshiftdb.org.

For example, the 13C-NMR spectra for butane is available at https://www.nmrshiftdb.org/NmrshiftdbServlet/nmrshiftdbaction/searchorpredict/smiles/CCCC/spectrumtype/13C.

Structured Text Files – JSON

JSON files are another common format for data storage and transmission. In R, the jsonlitepackage offers functions to import and convert this files to list or data frames.

if(!require("jsonlite")) {
  install.packages("jsonlite",
                   repos="https://cloud.r-project.org/")
  library("jsonlite")
}

jsBz <- fromJSON(
  "https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/241/JSON")
jsBz_NI <- jsBz$Record$Section[
  jsBz$Record$Section$TOCHeading=="Names and Identifiers",]
jsBz_NI_OI <- jsBz_NI$Section[[1]][
  jsBz_NI$Section[[1]]$TOCHeading=="Other Identifiers",]
jsBz_NI_OI$Section[[1]]$TOCHeading
##  [1] "CAS"                            "Related CAS"                   
##  [3] "Deprecated CAS"                 "European Community (EC) Number"
##  [5] "ICSC Number"                    "NSC Number"                    
##  [7] "RTECS Number"                   "UN Number"                     
##  [9] "UNII"                           "DSSTox Substance ID"           
## [11] "Wikipedia"                      "NCI Thesaurus Code"
CAS <- jsBz_NI_OI$Section[[1]][jsBz_NI_OI$Section[[1]]$TOCHeading=="CAS",]
table(unlist(CAS$Information[[1]]$Value$StringWithMarkup))
## 
## 26181-88-4    71-43-2 
##          1         13

Computational Programs – Excel

Let’s use the wine data set again.

str(wine)
## 'data.frame':    21 obs. of  31 variables:
##  $ Label                        : Factor w/ 3 levels "Saumur","Bourgueuil",..: 1 1 2 3 1 2 2 1 3 1 ...
##  $ Soil                         : Factor w/ 4 levels "Reference","Env1",..: 2 2 2 3 1 1 1 2 2 3 ...
##  $ Odor.Intensity.before.shaking: num  3.07 2.96 2.86 2.81 3.61 ...
##  $ Aroma.quality.before.shaking : num  3 2.82 2.93 2.59 3.43 ...
##  $ Fruity.before.shaking        : num  2.71 2.38 2.56 2.42 3.15 ...
##  $ Flower.before.shaking        : num  2.28 2.28 1.96 1.91 2.15 ...
##  $ Spice.before.shaking         : num  1.96 1.68 2.08 2.16 2.04 ...
##  $ Visual.intensity             : num  4.32 3.22 3.54 2.89 4.39 ...
##  $ Nuance                       : num  4 3 3.39 2.79 4.04 ...
##  $ Surface.feeling              : num  3.27 2.81 3 2.54 3.38 ...
##  $ Odor.Intensity               : num  3.41 3.37 3.25 3.16 3.54 ...
##  $ Quality.of.odour             : num  3.31 3 2.93 2.88 3.36 ...
##  $ Fruity                       : num  2.88 2.56 2.77 2.39 3.16 ...
##  $ Flower                       : num  2.32 2.44 2.19 2.08 2.23 ...
##  $ Spice                        : num  1.84 1.74 2.25 2.17 2.15 ...
##  $ Plante                       : num  2 2 1.75 2.3 1.76 ...
##  $ Phenolic                     : num  1.65 1.38 1.25 1.48 1.6 ...
##  $ Aroma.intensity              : num  3.26 2.96 3.08 2.54 3.62 ...
##  $ Aroma.persistency            : num  2.96 2.81 2.8 2.58 3.3 ...
##  $ Aroma.quality                : num  3.2 2.93 3.08 2.48 3.46 ...
##  $ Attack.intensity             : num  2.96 3.04 3.22 2.7 3.46 ...
##  $ Acidity                      : num  2.11 2.11 2.18 3.18 2.57 ...
##  $ Astringency                  : num  2.43 2.18 2.25 2.18 2.54 ...
##  $ Alcohol                      : num  2.5 2.65 2.64 2.5 2.79 ...
##  $ Balance                      : num  3.25 2.93 3.32 2.33 3.46 ...
##  $ Smooth                       : num  2.73 2.5 2.68 1.68 3.04 ...
##  $ Bitterness                   : num  1.93 1.93 2 1.96 2.07 ...
##  $ Intensity                    : num  2.86 2.89 3.07 2.46 3.64 ...
##  $ Harmony                      : num  3.14 2.96 3.14 2.04 3.64 ...
##  $ Overall.quality              : num  3.39 3.21 3.54 2.46 3.74 ...
##  $ Typical                      : num  3.25 3.04 3.18 2.25 3.44 ...

To write a data.frame to Excel, we can use the write_xlsx from the writexl package. To read an Excel file, read_xlsx from the readxl package can be used.

if(!require("writexl")) {
  install.packages("writexl",
                   repos="https://cloud.r-project.org/")
  library("writexl")
}

if(!require("readxl")) {
  install.packages("readxl",
                   repos="https://cloud.r-project.org/")
  library("readxl")
}
write_xlsx(list(wineSheet = wine),
                  path = "wine.xlsx")
wineExcel <- readxl::read_xlsx(path = "wine.xlsx")

Computational Programs – Matlab

To open and write Matlab data files, we can use the R.matlab package. Relevant functions are called readMat and writeMat.

if(!require("R.matlab")) {
  install.packages("R.matlab",
                   repos="https://cloud.r-project.org/")
  library("R.matlab")
}
writeMat("wine.mat", wine = wine)
wineMat <-readMat("wine.mat")
wineMat <- data.frame(as.character(unlist(wineMat[[1]][[1]])),
                      as.character(unlist(wineMat[[1]][2])),
                      as.data.frame(wineMat[[1]][3:31]))

Reading a Matlab data file returns a list that will need ad hoc manipulation.

RData Files

To save and read data in the R data format, we use save and load. These allow storing and restoring any set of variables of the environment. When loaded, variables are recovered with the same names they had on saving.

save(wine, file="wine.rda")
print(load("wine.rda"))          # print shows the name of the loaded objects
wineR <- get(load("wine.rda"))   # get allows storing the loaded information 
                                 # with a different name

Chemical Data Formats

Besides the general data file formats already discussed, there are several formats used specifically to store chemical information. Some significant ones are

More information on chemical data files can be found at https://en.wikipedia.org/wiki/Chemical_file_format, https://en.wikipedia.org/wiki/Mass_spectrometry_data_format, http://unichrom.com/chrom/uc-ffe.shtml or https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3518119/.

Chemical Data Formats – Chemical Table File

Chemical Table Files are a set of text-based file types designed to store molecular information, e.g. positions of the atoms and connection tables. Some of this format allow storing additional information such as conformers, additional properties and identifiers, or spectra.

Currently the simplest way to read a MOL or SDF file in R is the read.SDFset function in the ChemmineR Bioconductor package.

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager",
                   repos="https://cloud.r-project.org/")

if (!require("ChemmineR", quietly = TRUE)) {
  BiocManager::install("ChemmineR", update=FALSE)
  library("ChemmineR")
}
sdf1 <- read.SDFset(
  "https://cactus.nci.nih.gov/chemical/structure/CCCC/file?format=sdf")
draw_sdf(sdf1[[1]], filename=NULL)

header(sdf1)
## $CMP1
##                                    Molecule_Name 
##                                          "C4H10" 
##                                           Source 
## "APtclcactv01042217252D 0   0.00000     0.00000" 
##                                          Comment 
##                                              " " 
##                                      Counts_Line 
##        " 14 13  0  0  0  0  0  0  0  0999 V2000"
MW(sdf1[[1]])
##     CMP 
## 58.1222

Chemical Data Formats – JCAMP-DX

JCAMP-DX is an text-based open standard to store and distribute spectral data. Currently the best way to read these files in R is the readJDX function in the readJDX package.

if (!require("readJDX")) {
  install.packages("readJDX",
                   repos="https://cloud.r-project.org/")
  library("readJDX")
}

Acetone IR

jdx1 <- readJDX ("../data/67-64-1-IR.jdx")
plot(jdx1$Acetone$x,jdx1$Acetone$y, type="l",
     xlab="wave number", ylab="T")

JCAMP-DX downloaded from https://webbook.nist.gov/cgi/cbook.cgi?ID=C67641, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, Eds. P.J. Linstrom and W.G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 20899, https://doi.org/10.18434/T4D303, (on January 5, 2022).

Chemical Data Formats – NIST MSP

MSP is a text-based NIST-promoted format to store collections of spectra. It’s a common format in the metabolomics community.

Downloadable spectral databases can be found at

While a functional package to read MSP files doesn’t seem to be available, these files can be read with standard functions for text files.

    msp1 <- readLines("../data/MSMS-Neg-MassBankEU.msp")
    msp1 <- paste(msp1, collapse="\n")
    msp1 <- unlist(strsplit(msp1, "\n\n"))

    msp1_1 <- unlist(strsplit(msp1[1],"\n"))
    spectrum <- msp1_1[1:length(msp1_1) > which(substr(msp1_1,1,10)=="Num Peaks:")]
    spectrum <- read.table(text=spectrum, sep="\t")
    colnames(spectrum) <- c("m_z","int")
    spectrum
##        m_z  int
## 1  96.9602   70
## 2 241.0031 1000
    spectrum <- data.frame(m_z=rep(spectrum$m_z,each=3),
                          int=rep(spectrum$int,each=3),
                          i=1:3)
    spectrum$int[spectrum$i!=2] <- 0

    plot(x=spectrum$m_z, y=spectrum$int, xlim=c(40,300),
         type="l", xlab="m/z", ylab="")

::: {.bibref} MSMS-Neg-MassBankEU.msp downloaded from http://prime.psc.riken.jp/compms/msdial/main.html#MSP

Manipulación avanzada de tablas de datos

Tabla de datos o data frame

Com se ha viso anteriormente, el data frame es el tipo de dato más usado para almacenar tablas de datos.

A menudo, para poder realizar gráficos y/o aplicar distintos procedimientos estadísticos, es necesario manipular la tabla de datos. A esto dedicaremos este apartado.

Entre las operaciones habituales que haremos sobre una tabla de datos, hay

  • renombrar filas o columnas,
  • añadir columnas o filas,
  • segmentar,
  • cambiar el formato de un conjunto de datos,
  • eliminar filas o columnas,
  • crear tablas de datos de resumen, y
  • unir tablas

Terminaremos este bloque introduciendo los funciones básicas del paquete dplyr que facilita la realización de algunas de estas operaciones.

Partimos de una tabla de datos sintética…

df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
df
##   X1.5 letters.1.5. c.rep..a...3...rep..b...2..
## 1    1            a                           a
## 2    2            b                           a
## 3    3            c                           a
## 4    4            d                           b
## 5    5            e                           b

Renombrar filas o columnas

colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Añadir columnas o filas

df2 <- cbind(df, rnorm(5)) # añadir un vector al data frame
df2$var5 <- 5:1 # assignando valores a una nueva variable
df2
##            var1 var2 var3   rnorm(5) var5
## subject001    1    a    a -0.5952806    5
## subject002    2    b    a  1.6741805    4
## subject003    3    c    a -0.9135890    3
## subject004    4    d    b -1.5143077    2
## subject005    5    e    b  1.7074426    1
df2 <- rbind(df, list(6, "e", "b"))
df2
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b
## 1             6    e    b

Segmentar

Existen tres formas básicas de seleccionar filas o columnas de una tablas de datos

  • mediante índices numéricos,
  • usando los nombres de filas y columnas, y
  • mediante vectores lógicos

Segmentar – mediante índices

df[1:3,]
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
df[,c(1,3)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject003    3    a
## subject004    4    b
## subject005    5    b
df[-3,-2]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject004    4    b
## subject005    5    b

Segmentar – usando nombres

df[,"var2"]
## [1] "a" "b" "c" "d" "e"
df$var3
## [1] "a" "a" "a" "b" "b"
df[,c("var2","var3")]
##            var2 var3
## subject001    a    a
## subject002    b    a
## subject003    c    a
## subject004    d    b
## subject005    e    b

Segmentar – mediante vectores lógicos

df[c(T,T,F,T,F), c(T,F,T)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject004    4    b
df[df[,1] == 3 | df[,3] == "b",]
##            var1 var2 var3
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Cambiar el formato de un conjunto de datos

Un mismo conjunto de datos puede representar en una tabla de datos de acuerdo con distintas oragnizaciones o formatos.

Se denomina formato ordenado, tidy, aquella organización en la que las filas representan individuos, las columnas, variables de medidas, y las intersecciones los valores de dichas medidas.

Cuando existen más de una variable que corresponde al mismo tipo de medida, entonces los datos pueden organizarse de dos formas

  • usando una columna distinta para cada variable, formato ancho y tidy, o
  • usando una columna para los valores y otra para indicar cuál es la variable representada en esta fila.

El formato más adecuado dependerá de los análisis y visualizaciones que quieran hacerse.

Partiremos de un subconjunto de flights…

if(!require("nycflights13")) {
  install.packages("nycflights13")
  library("nycflights13")
}
fl_ny2ws_W <- flights[flights$dest %in% c("IAD","BWI"),
                    c("origin","dest","carrier","arr_delay","dep_delay")]
fl_ny2ws_W <- cbind(key = 1:nrow(fl_ny2ws_W), fl_ny2ws_W)
head(fl_ny2ws_W)
##   key origin dest carrier arr_delay dep_delay
## 1   1    LGA  IAD      EV       -14        -3
## 2   2    LGA  BWI      WN       -19        -1
## 3   3    EWR  IAD      EV        12        24
## 4   4    JFK  BWI      MQ       851       853
## 5   5    JFK  IAD      B6        14        15
## 6   6    LGA  IAD      EV         5        10

Cambiar el formato de un conjunto de datos – ancho a largo

fl_ny2ws_L <- rbind(
  cbind(edge = rep("origin", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
        airport = fl_ny2ws_W[,2], delay = fl_ny2ws_W[,6]),
  cbind(edge = rep("dest", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
        airport = fl_ny2ws_W[,3], delay = fl_ny2ws_W[,5]))
colnames(fl_ny2ws_L) <- c("edge","key","carrier","airport","delay")
head(fl_ny2ws_L)
##     edge key carrier airport delay
## 1 origin   1      EV     LGA    -3
## 2 origin   2      WN     LGA    -1
## 3 origin   3      EV     EWR    24
## 4 origin   4      MQ     JFK   853
## 5 origin   5      B6     JFK    15
## 6 origin   6      EV     LGA    10
tail(fl_ny2ws_L)
##       edge  key carrier airport delay
## 14957 dest 7476      B6     IAD    -4
## 14958 dest 7477      MQ     BWI    -9
## 14959 dest 7478      EV     IAD   -16
## 14960 dest 7479      EV     IAD   -19
## 14961 dest 7480      9E     IAD   -40
## 14962 dest 7481      9E     BWI     2

Cambiar el formato de un conjunto de datos – largo a ancho

fl_ny2ws_W2p1 <- fl_ny2ws_L[fl_ny2ws_L$edge=="origin",]
fl_ny2ws_W2p2 <- fl_ny2ws_L[fl_ny2ws_L$edge=="dest",]
fl_ny2ws_W2p1 <- fl_ny2ws_W2p1[order(fl_ny2ws_W2p1$key),-1]
fl_ny2ws_W2p2 <- fl_ny2ws_W2p2[order(fl_ny2ws_W2p2$key),-c(1,2)]
fl_ny2ws_W2 <- cbind(fl_ny2ws_W2p1,fl_ny2ws_W2p2[,-1])
colnames(fl_ny2ws_W2) <- c("key", "carrier", "origin", "dep_delay",
                           "dest","arr_delay")
head(fl_ny2ws_W2)
##   key carrier origin dep_delay dest arr_delay
## 1   1      EV    LGA        -3  IAD       -14
## 2   2      WN    LGA        -1  BWI       -19
## 3   3      EV    EWR        24  IAD        12
## 4   4      MQ    JFK       853  BWI       851
## 5   5      B6    JFK        15  IAD        14
## 6   6      EV    LGA        10  IAD         5
tail(fl_ny2ws_W2)
##       key carrier origin dep_delay dest arr_delay
## 7476 7476      B6    JFK        12  IAD        -4
## 7477 7477      MQ    JFK        -5  BWI        -9
## 7478 7478      EV    LGA        -5  IAD       -16
## 7479 7479      EV    JFK        -1  IAD       -19
## 7480 7480      9E    JFK        -7  IAD       -40
## 7481 7481      9E    JFK        38  BWI         2

Eliminar filas o columnas

La forma más habitual de eliminar filas o columnas es segmentando la tabla de datos. Sin embargo, una columna también puede eliminarse asignando la misma a NULL.

df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df <- df[-2,]
df$var2 <- NULL
df
##            var1 var3
## subject001    1    a
## subject003    3    a
## subject004    4    b
## subject005    5    b

Si lo que se desea es eliminar una variable del entorno de trabajo entonces se usa la función rm.

rm(df)

Creación de resúmenes a partir de datos en formato ancho

sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))

sum_fl_ny2ws$mean_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                                 mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                      function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
                              sd,na.rm=TRUE)
sum_fl_ny2ws
##     edge mean_delay r_delay  s_delay
## 1 origin   16.84267     885 47.85548
## 2   dest   13.11556     904 51.52259

Creación de una tabla de resumen a partir de datos en formato largo

sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))
sum_fl_ny2ws$mean_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
                              mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
                           function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,sd,na.rm=TRUE)
sum_fl_ny2ws
##     edge mean_delay r_delay  s_delay
## 1 origin   13.11556     904 51.52259
## 2   dest   16.84267     885 47.85548

dplyr

Estas mismas operaciones que se han comentado más arriba, se pueden llevar a cabo también a partir de un sistema coherente de métodos creado como una gramática para la manipulación de datos, el paquete dplyr –que forma parte de tidyverse–.

http://dplyr.tidyverse.org/

En dplyr las operaciones básicas, se realizan mediante cuatro métodos:

  • select: segmentar columnas,
  • filter: segmentar filas,
  • mutate: crear nuevas columnas, y
  • arrange: ordenar.

Estas se encadenan usando el operador pipe, %>%.

df <- flights %>% dplyr::select(origin, dest, arr_delay) %>% 
  filter(origin == "LGA" & (dest == "IAD" | dest == "BWI")) %>%
  mutate(arr_delay_h=arr_delay/60) %>% 
  arrange(-arr_delay_h)
df
## # A tibble: 1,818 x 4
##    origin dest  arr_delay arr_delay_h
##    <chr>  <chr>     <dbl>       <dbl>
##  1 LGA    IAD         398        6.63
##  2 LGA    IAD         326        5.43
##  3 LGA    IAD         312        5.2 
##  4 LGA    IAD         306        5.1 
##  5 LGA    IAD         292        4.87
##  6 LGA    IAD         281        4.68
##  7 LGA    IAD         280        4.67
##  8 LGA    IAD         265        4.42
##  9 LGA    IAD         262        4.37
## 10 LGA    IAD         253        4.22
## # ... with 1,808 more rows

Para la creación de resúmenes a partir de tablas en formato largo, es muy útil y cómoda la combinación group_by y summarize.

df %>% group_by(dest) %>% summarise(mean_delay = mean(arr_delay, na.rm=TRUE))
## # A tibble: 2 x 2
##   dest  mean_delay
##   <chr>      <dbl>
## 1 BWI        -7.87
## 2 IAD        13.7

Por último, las conversiones entre formatos también pueden hacerse a partir de las funciones spread y gather del paquete tidyr –incluido en tidyverse y que se integra de forma natural en la gramática propuestas por dplyr–.

Advanced Graphics in R

Gramática de gráficos (GoG)

La gramática de gráficos es una aproximación teórica al estudio de los componentes de un gráfico. De acuerdo con este análisis, un gráfico se puede construir mediante la especificación de un conjunto de capas y componentes que definen los datos, la asociación de los mismos a aspectos del gráfico, la especificación de la relación entre los valores de las variables de los datos con las del gráfico, la estructura geométrica del gráfico…

Wilkinson, L. (2006). The grammar of graphics. Springer Science & Business Media.

ggplot2

ggplot2 es un paquete de R que implementa de la gramática de gráficos.

Wickham, H. (2010). A layered grammar of graphics. Journal of Computational and Graphical Statistics, 19(1), 3-28.

Referencias:

ggplot2 forma parte del paquete tidyverse (aunque también puede instalarse y cargarse autónomamente).

if(!require("tidyverse")) {
  install.packages("tidyverse", repos="https://cloud.r-project.org/",
         quiet=TRUE, type="binary")
  library("tidyverse")
}

ggplot2 – capas y elementos del gráfico

En ggplot2 cada elemento gráfico que representa un conjunto de datos constituye una capa. Una o más capas constituyen un gráfico.

Cada capa queda definida mediante la especificación de sus elementos. Los principales son

  • datos,
  • mapeado estético, y
  • geometrías

En ggplot2, los gráficos constituyen un objeto de R y se construyen de forma aditiva.

Por ejemplo,

grafico <- ggplot(data = anscombe,
        mapping = aes(x = x1, y = y1))  # Datos y mapeado estético 
grafico <- grafico + geom_point()       # Geometría

grafico

ggplot2 – datos

En ggplot2, el elemento data (datos) se introduce como primer argumento de la función ggplot. Debe corresponder a una tabla de datos o un tipo de datos convertible a tabla de datos.

grafico <- ggplot(data = anscombe,

ggplot2 – mapeado estético

El mapping (mapeado estético) corresponde al establecimiento de relaciones entre variables de los datos y variables del gráfico. Es el segundo argumento de la función ggploty debe crearse con al función de apoyo aes.

        mapping = aes(x = x1, y = y1))

Para variables cuantitativas, los mapeados más comunes corresponden a

  • posiciones: x, y…
  • tamaño: size
  • color: color, fill

Para variable cualitativas, los mapeados más frecuentes son

  • posiciones: x, y…
  • color: color, fill
  • forma: shape

ggplot2 – geometrías

Las geometrías (geom_) indican la forma que debe tener el gráfico, es decir, cómo se articulan las variables del gráfico. Se añaden al gráfico sumándose al objeto creado por ggplot.

grafico <- grafico + geom_point()

Son geometrías de uso común

  • geom_point
  • geom_line, geom_vline, geom_hline
  • geom_bar
  • geom_histogram
  • geom_boxplot

Un resumen de las geometrías y su relación con las variables del gráfico que reconoce cada una de ellas figura en https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf.

Gráficos en ggplot2

Veremos cómo crear en ggplot2 los gráficos más habituales, añadiendo algunas consideraciones para aquellos casos donde los gráficos realizados con base tienen prestaciones insuficientes.

  • gráficos de dispersión,
  • histogramas,
  • diagramas de barras, y
  • diagramas de caja (boxplot).

De forma general, los gráficos en ggplot2 se construyen a partir de tablas de datos (data frames), de los cuales se seleccionan las variables a representar.

Usaremos 1000 datos del conjunto de datos diamonds para crear los distintos ejemplos.

diaM <- diamonds[sample(1:nrow(diamonds),1000),]
str(diaM)
## tibble [1,000 x 10] (S3: tbl_df/tbl/data.frame)
##  $ carat  : num [1:1000] 0.55 2.07 1.29 0.44 0.26 0.94 0.31 0.3 0.51 0.41 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 4 3 5 5 3 4 2 5 5 4 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 1 6 7 4 3 3 4 4 2 2 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 5 3 7 5 7 2 6 4 5 3 ...
##  $ depth  : num [1:1000] 62 63.1 61.2 61.9 60.6 62.3 63.3 62.5 62.8 62.8 ...
##  $ table  : num [1:1000] 58 59 56 54 56 59 57 53 57 58 ...
##  $ price  : int [1:1000] 1975 15751 6918 915 657 3095 707 526 2075 904 ...
##  $ x      : num [1:1000] 5.26 8.11 7.01 4.89 4.12 6.32 4.28 4.27 5.1 4.77 ...
##  $ y      : num [1:1000] 5.28 8.02 7.05 4.93 4.16 6.23 4.31 4.3 5.07 4.72 ...
##  $ z      : num [1:1000] 3.27 5.09 4.3 3.04 2.51 3.91 2.72 2.68 3.2 2.98 ...

Gráficos en ggplot2 – gráfico de dispersión

ggplot(diaM, aes(x=carat,y=price)) + geom_point()

Añadiendo una tercera variable (cut) y modificando algunos aspectos de formato…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3)

Añadiendo líneas de tendencia…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3) +
  geom_smooth(method="lm",se=FALSE)
## `geom_smooth()` using formula 'y ~ x'

Gráficos en ggplot2 – histograma

ggplot(diaM, aes(x=price)) + geom_histogram(binwidth=1000)

Y en función del corte…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

En frecuencias relativas (por grupo)…

ggplot(diaM, aes(x=price,y=..density..,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

Quizás funcione mejor un gráfico de densidades…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_density(alpha=.3)

Gráficos en ggplot2 – diagrama de barras

ggplot(diaM, aes(x=clarity)) + geom_bar()

En función de la claridad…

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar()

Para comparar entre frecuencias absolutas, funcionan mejor las barras separadas.

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="dodge")

Para comparar entre frecuencias relativas acumuladas, son mejores las barras apiladas en frecuencia relativa (para cada clase).

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="fill")

Los diagramas de barras también pueden crearse a partir de tablas de datos agrupados. En este caso, debe indicarse qué variable es la y e incluir stat="identity" en el geom_bar.

Gráficos en ggplot2 – diagrama de caja

ggplot(diaM, aes(x=1, y=price)) + geom_boxplot()

NOTA: Para el geom_boxplot se requieren dos variables. Si no hay variable independiente se puede incluir x=1.

Y en función del corte…

ggplot(diaM, aes(x=cut, y=price)) + geom_boxplot()

El gráfico se puede mejorar mostrando todos los puntos, con una posición aleatorizada.

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(shape = 21, alpha=.5,height=0,width=.2)

O incluyendo un violin plot y un punto para la media…

medias <- diaM %>% group_by(cut) %>%
  summarise(price=mean(price))

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.1) +
  geom_point(data=medias,shape=3)

ggplot2 – elementos adicionales

Además de los elementos ya presentados (datos, mapeado estético y geometrías), otros elementos de ggplot2 permiten controlar aspectos adicionales del gráfico, por ejemplo

  • scale_...: controlan los aspectos relativos a la presentación de las variables del gráfico,
  • coord_...: establecen el sistema de coordenadas usado en la geometría,
  • labs: establece los títulos del gráfico,
  • theme, theme_...: controlan la part del gráfico que no corresponde a datos (non-data ink), y
  • facet_: permite la creació de secuencia de gráficos en función de una o dos variables

Un ejemplo para terminar…

ggplot(diaM, aes(x=carat, y = price, shape = cut, col = clarity)) +
  geom_point(alpha=.6) +
  scale_x_continuous(breaks=1:3) +
  scale_y_continuous(trans="log10") +
  scale_color_brewer(palette="Spectral")+
  facet_grid(cut ~ clarity) + 
  theme_bw() + 
  theme(legend.position = "none",text = element_text(size=10))
## Warning: Using shapes for an ordinal variable is not advised